## To do: test these functions!!

## To run SRE model over a range of range-edge gradients, dispersal costs, and Repro rates

source(file="/home1/30/jc227089/evo-dispersal/SRE/Stable range functions.R", echo=F)
setwd("/home1/30/jc227089/SRE/gvar/ran_outs")
load("/home/30/jc227089/SRE/equil/ran_outs/SRE-equil_summary.RData")
### Functions ###

# to collect data
collect<-function(popmatrix, edge.X, last.occ, spX, spY, gen, ...){ #edge x is last sustainable population column
	if (length(popmatrix[,1])==0) return(NULL)	
	popsize<-tapply(popmatrix[,"P"], popmatrix[,"X"], length) #popsize over columns
	mean.disp<-tapply(popmatrix[,"P"], popmatrix[,"X"], mean) #mean trait value over columns
	occpd<-table(factor(popmatrix[,"Y"], levels=1:spY), factor(popmatrix[,"X"], levels=1:spX))>0
	occpd<-occpd[,(edge.X+1):length(occpd[1,])] # occupied cells to R of edge.X
	n.occpd<-sum(occpd) # number occupied to R of edge.X
	tover<-occpd!=last.occ
	tover<-sum(tover) # turnover of sites R of edge.X
	gen<-rep(gen, length(popsize))
	out<-list(c(..., n.occpd, tover), cbind(gen, popsize, mean.disp))
	out
}

gvar<-function(n=20*30*20, ngens.init=250, ngens=500, maxX=30, spY=20, end.samp=20, steps, K=40, lambda, d.cost, mutn=0.05, plot=T, ep.array){
	pop<-init.inds(n, maxX, spY)
	hab<-init.hab(maxX, spY, 0)
	nbrs<-neighbours.init(maxX, spY) #reveal all neighbours
	edge<-edge.finder(lambda, d.cost, ep.array, maxX, steps) ##lookup edge X based on lambda and d.cost
	data<-list(par=c(), space=c()) # to take last 50 gens
	for (i in 1:ngens.init){
		pop<-repro.disp(pop, hab, maxX, spY, K, lambda, d.cost, nbrs, mutn)
		if (plot==T) plotter(pop, maxX, spY, K)
	}
	hab<-init.hab(maxX, spY, steps) #drop habitat gradient over the population
	for (i in (ngens.init+1):ngens){
		pop<-repro.disp(pop, hab, maxX, spY, K, lambda, d.cost, nbrs, mutn)
		if (plot==T) plotter(pop, maxX, spY, K)
		if (ngens-i==(end.samp+1)) {
			occpd<-table(factor(pop[,"Y"], levels=1:spY), factor(pop[,"X"], levels=1:maxX))>0
			occpd<-occpd[,(edge+1):length(occpd[1,])] # occupied cells to R of edge.X
		}
		if (ngens-i<=end.samp) {
			temp<-collect(pop, edge, occpd, maxX, spY, i, steps, lambda, d.cost, edge)
			data[[1]]<-rbind(data[[1]], temp[[1]])
			data[[2]]<-rbind(data[[2]], temp[[2]])
			occpd<-table(factor(pop[,"Y"], levels=1:spY), factor(pop[,"X"], levels=1:maxX))>0
			occpd<-occpd[,(edge+1):length(occpd[1,])] # occupied cells to R of edge.X
		}
	}
	colnames(data[[1]])<-c("steps", "lambda", "d.cost", "edge", "n.occpd", "tover")
	xloc<-as.numeric(row.names(data[[2]]))
	data[[2]]<-cbind(xloc, data[[2]])
	data
}

### initialise ###
args=(commandArgs(TRUE))

#evaluate the arguments
# input arguments will be number of replicates runs (rr), file.ID, lambda (lv), steps (stps) and dispersal cost (dc)
for(i in 1:length(args)) {
	 eval(parse(text=args[[i]]))
}

pars<-c(lv, stps, dc)
sim.result<-c()
for (rep in 1:rr){
	out<-gvar(n=20*30*20, ngens.init=250, ngens=500, maxX=30, spY=20, end.samp=20, steps=stps, K=40, lambda=lv, d.cost=dc, mutn=0.05, plot=T, ep.array=sumry)
	out<-c(list(pars), out)
	sim.result<-cbind(sim.result, out)
}
save(sim.result, file=paste("gvarout", file.ID, "RData", sep=""))